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We develop a first-principles scheme to calculate adiabatic and non-adiabatic phonon frequencies 
in the full Brillouin zone. The method relies on the variational properties of a force-constants 
functional with respect to the first-order perturbation of the electronic charge density and on the 
localization of the deformation potential in the Wannier function basis. This allows for calculation 
of phonon dispersion curves free from convergence issues related to Brillouin zone sampling. In 
addition our approach justify the use of the static screened potential in the calculation of the 
phonon linewidth due to decay in electron-hole pairs. We apply the method to the calculation of 
the phonon dispersion and electron-phonon coupling in MgB2 and CaCe. In both compounds we 
demonstrate the occurrence of several Kohn anomalies, absent in previous calculations, that are 
manifest only after careful electron and phonon momentum integration. In MgB2, the presence 
of Kohn anomalies on the E25 branches improves the agreement with measured phonon spectra 
and affects the position of the main peak in the Eliashberg function. In CaCe we show that the 
non-adiabatic effects on in-plane carbon vibrations are not localized at zone center but are sizable 
throughout the full Brillouin zone. Our method opens new perspectives in large-scale first-principles 
calculations of dynamical properties and electron-phonon interaction. 

PACS numbers: 63.20.dk, 71.15.-m, 63.20.kd 
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I. INTRODUCTION 

Electron-phonon (EP) interaction is responsible of 
many important phenomena in solids. As an example, 
the temperature behavior of the electron relaxation time 
in metals is to a great extent due to the scattering be- 
tween carriers and atomic vibrations^ such that finite 
temperature transport is largely ruled by the EP inter- 
action. Similarly, the high temperature heat capacity in 
metals is enhanced by the increased electronic mass due 
to the interaction with lattice vibrations^. In metals, at 
low temperatures, EP coupling can generate a supercon- 
ducting state in which electrons move with no electrical 
resistance'^. It also can increase the effective mass of the 
carriers so much that the system is driven from a metallic 
to an insulating state, as it happens in case of polaronic or 
Peierls instabilities*. Finally, the electron-phonon scat- 
tering is often the largest source of phonon damping in 
phonon-mediated superconductors'^ . 

First-principles theoretical determination of the 
electron-phonon coupling strength in solids requires 
the calculation of the electronic structure, the vibra- 
tional properties, and the electron-phonon coupling 
matrix elements. In state-of-the-art electronic struc- 
ture calculations these quantities are obtained using 
the adiabatic Born-Oppenheimer approximation^ and 
density-functional theory (DFT) in the linear-response 
approach ^"^ . 

More specifically, within the Born-Oppenheimer ap- 
proximation, the determination of phonon frequencies, 
that are related to the real part of the phonon self- 
energy, requires the calculation, in a self- consistent man- 
ner, of the variation of the Kohn-Sham potential Vscf 
with respect to a static ionic displacement u^, namely 



— ^?^ ■ As the displacement of the ions is static, 

the obtained — ^^^ ■'^' is real. Conversely, the phonon 
linewidth, related to the imaginary part of the phonon 
self-energy, is obtained in a non self-consistent proce- 
dure using the electron charge density and the previously 
determined — ^^^^""^ . The advantage of using a non self- 
consistent procedure to study phonon linewidth is mainly 
related to the less expensive computational load with re- 
spect to a self-consistent one. In addition, as recently 
demonstrated, interpolation schemes^'^'^^ of the electron- 
phonon matrix elements can be used to calculate the 
imaginary part of the phonon self-energy on ultra-dense 
k-point grid. 

It is however unclear to what extent this procedure of 
calculating self-consistently the real part of the phonon 
self-energy and non self-consistently the imaginary part 
is actually correct. Indeed, the proper way of treating 
phonons should be to consider a monochromatic time- 
dependent displacement (at the phonon frequency w) and 
perform a time-dependent self-consistent linear-response 
scheme. The resulting variation of the self-consistent 
potential — ^'^/u would then be a complex quantity. 
The real part of the resulting phonon self-energy would 
then determine the phonon frequencies while its imagi- 
nary part would lead to the phonon linewidth. In this 
way, non-adiabatic (in the sense of Ref. 12) dynamical 
phonon frequencies could be accessed. Although feasible 
in principle this procedure would require a full rewriting 
of the linear response code including time dependence 
and it would also be more expensive then a standard 
static linear-response calculation. Moreover, in the pres- 
ence of Kohn anomalies and long-range force-constants, 
where extremely accurate k-point sampling of the Fermi 



surface is needed to converge the phonon self-energy, the 
calculation would be unfeasible. 

Thus, it would be desirable to have a non -self consis- 
tent linear-response formulation to obtain both the real 
and imaginary phonon self-energy, both within the adia- 
batic/static and non-adiabatic/dynamic approximation. 

In this work we develop a scheme to calculate non self- 
consistcntly both the real and the imaginary part of the 
phonon self-energy by using a functional that is varia- 
tional with respect to the variation of the self-consistent 
charge density. Our method opens the way to calculate 
adiabatic/static and non-adiabatic/dynamic phonon fre- 
quencies using an ultra-dense sampling of the electron 
and phonon wave-vector in a non self-consistent way, 
starting from the static self-consistent variation of the 

time independent Kohn-Sham potential ( — ^^u ) ' '-'^" 
tained with a coarse sampling. In this way, the main com- 
putational load related to phonon frequencies calculation 
is drastically reduced. The efficiency of the method is fur- 
ther enhanced by using the interpolation of the electron- 
phonon matrix elementS"'^^ based on Wannier functions^^. 

The method is applied to study the dynamical and 
superconducting properties of MgB2 and CaCg, two of 
the most studied superconducting materials in the last 
few years. These applications are meaningful and com- 
putationally challenging. In fact, although many exper- 
imental and theoretical studies appeared, there are still 
important open issues and debated results. 

For example, conflicting results for the calculate MgB2 
electron-phonon coupling (A) are present in the litera- 
ture mainly due to difficulties in Brillouin Zone sampling. 
This, apart its fundamental importance, prevents a full 
understanding of the normal and superconducting prop- 
erties of this material. 

In the case of CaCg, very large non-adiabatic effects 
were predicted^^ and measured^^ at zone-center. It is 
unclear to what extend non-adiabatic effects are sizable 
far from F-point. In this last case, the adiabatic effects 
can be relevant for thermodynamic properties resulting 
from an average of the phonon frequencies over the Bril- 
louin zone. 

The paper is organized as follows. In Section II we 
derive the variational formulation of the force-constant 
matrix. This will be done starting from the formal def- 
inition of force constants in the case of a monochro- 
matic time-dependent ionic displacement (Sec. II A). In 
Sec. II B we discuss when non-adiabatic/dynamic effects 
should be expected in the phonon spectra. The linear 
response equations for the dynamical matrix are intro- 
duced (Sec. II C) in the general case and then specialized 
in the density functional formulation (Sec. II D). 

In Section II E, we develop the variational formulation 
of the force constants in the density functional linear- 
response theory and outline the computational frame- 
work of phonon and electron-phonon coupling (Sec. II F 
and II G). In Sec. II H we show how our approach justify 
the use of the static screened potential in the calculation 
of the phonon linewidth due to decay in electron-hole 



pairs. Section III will be devoted to the description of 
the implementation of the theory using the Wannier in- 
terpolation scheme. 

In Section IV we reports results on dynamical and su- 
perconducting properties of MgB2 (Sec. IV A) and CaCe 
(Scc.IVB). 

Section V summarizes our conclusions. 



II. THEORY 

A. Time dependent dynamical matrix and phonon 
damping. 

We consider a crystal with N atoms in the unit cell 
submitted to a time dependent perturbation, in which 
the position of an atom is identified by the vector 



ivr ^ R,^ 



U7(f), 



(1) 



where Rl is the position of the L-th unit cell in the 
Bravais lattice, r <, is the equilibrium position of the s-th 
atom in the unit cell, u/(i) indicates the deviation from 
equilibrium of the nuclear position and / = {L, s}. The 
force at time t acting on the J-th nucleus (J == {M, r}) 
due to the displacement u/(t') of the atom I-th at time 
t' is labeled Fj{t). The force constants matrix is defined 
as: 



Cij (Rl — Rm ',t — t) 



SFjjt) 
Suiit') 



(2) 



where we used the translational invariance of the crystal 
and make evident the dependence of Cjj on the lattice 
vector Rl — Rm (to lighten the notation we omit it in 
the following equations where no confusion may arise). 
The w-transform of the force-constants matrix is thus: 



(3) 



Cuiu) = / dte^'^'Cijit) 



While the force-constants matrix Cij{t) is a real quan- 
tity, its cj-transform Cij{uj) is not real and has both a 
real and imaginary part. The Fourier transform of the 
force-constant matrix is 



C,,(q,cu) = ^e-*'i^-CL,.Mr(c< 



(4) 



where, without loss of generality, we have chosen Rm ~ 
0. The Hermitian and anti-Hermitian combination of the 
force-constant matrix in momentum space are: 



£'sr(q,tj) 



1 



2^/MsMr 
1 



2iy/MsMr 



[CsM,Uj)+Crs{c{.i^)*] (5) 

[CsM,oo)-Crs{ci,u:)*] (6) 



where Mg is the mass of the s-th atom in the unit cell. 
These quantities are associated to the real and imaginary 



part of the dynamical matrix in coordinate space. If the 
imaginary part of the dynamical matrix is small with 
respect to its real part, namely 



\Asr{q,uj)\ « \Dsr{q,uj) 



then the self-consistent condition 



det \Dsr{ci,UJqu) - ^IJ = 



(7) 



(8) 



determines non-adiabatic/dynamic phonon frequencies 



w, 



j^ and phonon eigenvectors {eq^} and ly 



1,3N indicates the phonon branches. The adia- 
batic/static phonon frequencies and eigenvectors are ob- 
tained considering a static perturbation, thus diagonal- 



izmgDrs{(l, 



0). 



On the other hand, the imaginary part of the force- 
constants matrix determines the phonon-damping 



Iqu 






(9) 



and the phonon lincwidth. 



B. Adiabatic and non-adiabatic phonons 

In the previous section we have defined the non- 
adiabatic/dynamic phonon frequencies as the eigenval- 
ues of the Hermitian part of the time dependent dynam- 
ical matrix, and the adiabatic/static phonon frequen- 
cies as the eigenvalues of the static time-independent 
dynamical matrix. Solid-state text-books"'^'^^ and first- 
principles calculations of the phonon dispersion'''^'^^"^^, 
usually treat only the adiabatic/static case, since it is 
commonly assumed that the adiabatic phonon frequen- 
cies coincide with the non-adiabatic ones. 

In insulators, where the fundamental energy gap be- 
tween the electronic ground state and the first available 
excited state is much larger than the phonon energy, the 
adiabatic/static approximation is well justified. In met- 
als the situation is more complex. ^^'^"'^^ 

The crucial parameter in metals is the electron relax- 
ation time T. In absence of electron-defect, electron- 
electron and electron-phonon scattering, r is infinite. In 
a real metallic system, the presence of these scattering 
processes results in a finite relaxation time r. We can 
define three cases: (i) a clean limit, when the electron 
relaxation time t is much larger than the phonon pe- 
riod divided by 27r (ii) a dirty limit, when the electron 
relaxation time r is much smaller than the phonon pe- 
riod divided by 27r (iii) an intermediate regime, when the 
electron relaxation time r is comparable to the phonon 
period divided by 27r. 

It has been shown^^^^^'^^ that in the dirty limit the 
non-adiabatic phonon frequencies coincide with the adi- 
abatic ones also in metals. Instead, in the clean limit 
and in the intermediate regime, the adiabatic and non- 
adiabatic frequencies are, in general, different. 



The non-adiabatic calculations based on time- 
dependent DFT (described in the following sections) are 
performed in the perfect clean limit. Indeed, in our time- 
dependent calculations, the electron relaxation time r is 
infinite, since we use an instantaneous (real in the fre- 
quency space) exchange-correlation Kernel, and we do 
not consider a broadening of the electronic levels due to 
defects and electron-phonon scattering. 

Thus, our non-adiabatic/dynamic DFT frequencies 
should be used to reproduce the phonon frequencies mea- 
sured in metal in the clean-limit. Instead, the adia- 
batic/static DFT frequencies (those generally computed 
with DFT linear response codes) should be used to repro- 
duce the phonon frequencies measured in a metal in the 
dirty limit. Finally, to reproduce the phonon frequen- 
cies measured in a metal in the intermediate regime, one 
should explicitly include electron-lifetime effects in the 
linear response calculation, as it has been done for zone 
center phonon in Refs. 12,22. 

By considering a single band and by linearizing the 
electronic-band dispersion near the Fermi energy, it has 
been shown^" that the differences between the adiabatic 
and non-adiabatic phonons are largest at the center of 
the Brillouin zone (BZ) and vanish for q ^ Wq^/wp, 
where v-p is the Fermi velocity. Such differences, at the 
BZ center and in the clean limit, have been computed 
within DFT^^, and can be very sizable (up to 30% of the 
phonon frequency). However a detailed DFT study of 
non-adiabatic effect away from the BZ center (beyond a 
linearized one band approximation^^) is still missing. 



C. Time dependent linear response theory 

The force-constant matrix can be evaluated in the lin- 
ear response theory, considering that the atomic displace- 
ment induces a perturbation in the external potential act- 
ing on the electrons. 

The Hellmann-Feynmann theorem^^"^^ states that the 
force on atom J, Fj(t) can be evaluated in terms of the 
variation of the external potential: 



Fj(t) 



drn{r,t) 



SVc.t{r) 



(10) 



where n{r,t) is the electronic charge density and T4xt(r) 
is the external potential, namely: 



K=xt(r) 



E 



z 



1 



ZrZ. 



I^J 



,r-R/| 2 ^^ IRj-R/l 



(11) 



It is worthwhile to recall that T4xt(r) does not depend 
explicitly on time but only through the dependence on 
time of the phonon displacement u/ (i) . In linear response 
the force-constant matrix Cjj{uj) is written as : 



Cij{lu) 



• Sn{r,uj) SVc^tjr) 
5ui 5uj 

/ "0 r ^ — 7 dr. 

J SujSuj 



(12) 



where rio(r) is the unperturbed charge density, n(r,uj) 
is the w-transform of the time-dependent charge den- 
sity n(r, t) and — piLEl and , '-'y^'"-' are real and time- 
independent quantities evaluated at the equilibrium po- 



sition of the nuclei. On the contrary 



(5u/ 



is complex. 



D. Time-dependent linear response in density 
functional theory 

In this section, we examine the linear-response calcu- 
lation of the real-space force-constants in the framework 
of density functional theory. The derivations are kept 
in real space as we want to keep the presence of imag- 
inary terms (damping) in the force-constants matrix as 
manifest. 

In order to compute n(r, cj), we assume a monochro- 
matic perturbation of the form 



u/(i) = ujH(e''"*+e-*"*) 



(13) 



where U7(a;) is real. In density functional theory the 
resulting monochromatic perturbing potential is the ex- 
ternal potential 14xt(r) in Eq. 11 . The derivative of the 
w-transform of the charge density with respect to a ionic 
displacement is written as^ 



n]{r,Lj,T) = 



dn{r,Lu) 



5ui 

{ip\^'j\SVscF{r,uj)/Suj\tpi,i) 



Nk{T) 



7Ai:,(r)^k',(r) (14) 

where k,k' label the crystal momentum, i,j are band in- 
dexes, the T is the electronic temperature (or broaden- 
ing) in the Fermi function /ki (T) , ry is an arbitrarily small 
positive real number, Nk{T) is the number of k-points 
needed to converge the sum in Eq. 14 at the electronic 
temperature T and the factor 2 (here and in the follow- 
ing) accounts for the spin-degeneracy. Finally, Vscf is 
the Kohn-Sham self-consistent potential. The quantity 
^^SCF(r, w)/5u7 is complex and given by 

—^-— = ^—+JK{v^r)n,{v,u:,T)dv. 

(15) 

where Kir, r') — , ,^^^,,^ is the kernel of the Hartree 
and Exchange and correlation functional, -Ehxc["']- As 
usual we have assumed the Hartree and exchange- 
correlation Kernel to be instantaneous (real in lo space). 
Substitution of Eq. 14 in Eq. 12 leads to: 

Nk(T) 

Cu{uj,T)^2 ^ {MT)-h',{T)) 

ki.k'j 

{tp\^'j\SVscF (r, w)/(5u7 IV'ki) (V'ki It^Kixt (r)/(5u,/ l^k'i) 



Ckj — Ek'j 



ir] 



no{r) 



SujSuj 



dr. 



(16) 



where, from now on, we explicitly indicate the depen- 
dence on the electronic temperature T. This expression 
of the force-constants matrix is normally used in stan- 
dard implementations of linear- response theory^. 

Further substitution of Eq. 15 in Eq. 16 gives the 
following alternative, but equivalent formulation for the 
force-constants matrix in linear response theory: 



Cij{cj,T)^2 



""y^ /k.(T)-M(r) 

ki.k'j 



Cki - ek'j" +i^ + iil 



^ ,, |(5VscF(r,w) (5Fscf(i-,w) 

X (V'k'il J |V'ki)(V'kz| 7 IV'k'j) 



dv no (r 



Sui 

(5'Kxt(r) 



5uj 



SujSuj 
n]j{r, to, T)K{v, v')n\{v' , u, T)dvdv' . 



(17) 



In Eq. 17 the term including the second derivative of 
the external potential is real whereas all the other terms 
are complex. The advantage of Eq. 17 is that it allows 
us to introduce a variational formulation of the force- 
constants matrix as it will be shown in the following. 



E. Variational formulation of the force-constants. 

We introduce the following force-constants functional 
Fij: 



Fu[p{v),p'{v),u,T] 



Nk{T) 



= 2E 



MT)~h',{T) 



X (Vk',|^%^+ /i^(r,r')p(r')dr'|7^k.) 



X (V'l 



i5u/ 

,5Kxt(r) 



Suj 

drno(r) 



+ / X(r,r')p'(r')rfr'|^k'j) 



5^V,^t{r) 



SujSuj 
- f f p{r)K{r,r')p'{r')drdr'. (18) 

With this definition the force-constants matrix reads as: 

Cuiio, T) = Fij[n]{v, c^, T), n^(r, uo, T),lu, T] (19) 

It is straightforward to note that the force-constant 
functional is quadratic with respect to p, namely: 

SFij[pir),p'{r),Lo,T] 



Sp{i 



= 0. 



p{r)=nj{r,ij,T),p'{r)=n^j{r,u:,T) 



(20) 



The same relation holds upon derivation with respect to 
P'- 



The consequence of Eq. 20 is that Cij{uj,T) is an 
extremal point and that a given error on n]{r,uj,T) or 
on nj{r,uj,T) affects the functional and the phonon fre- 
quencies only at second order. 

Our functional formulation is related to that used for 
dielectric tensors by Gonze et at in Ref. 26. The main 
difference is that, while the functional formulation of Ref. 
26 is variational with respect to the first-order perturba- 
tion of the wave-function, our formulation is variational 
with respect to the first-order perturbation of the elec- 
tronic charge density. 



F. Approximated force-constants functional. 

At this point we can exploit the results of the pre- 
ceding section in order to develop a method to accu- 
rately calculate phonon and electron-phonon properties. 
Within density functional theory the most precise force- 
constant matrix C'ij{uj,T) is obtained when T coincides 
with the physical temperature To {e.g. room tempera- 
ture) of the system. However, in a metal, this tempera- 
ture is prohibitively small and the self-consistent calcu- 
lation of Ci.j{uj,Tq) would be computationally unfeasi- 
ble as it requires a large number of k-points in order to 
properly converge the summation in Eq. 17. In addition, 
the self-consistent calculation of the imaginary part of 
n}(r,w,T) and of nj(r,w,T) at finite w 7^ requires to 
increase further the number of k-points with respect to 
standard static linear-response calculations of n}(r, 0, T). 

A remedy to these problems is to define an approxi- 
mated force constants matrix, Cij^ as: 



a Fourier transform to obtain: 



C/j(c^,ro) = f7j[n}(r,0,rph),n^(r,0,Tph),c<^,ro] 



(21) 



that is, (i) the frequency dependence on n}(r, cj,Tph) is 
neglected and its static limit (w ~ 0) is considered and 
(ii) n} and tVj are calculated at the temperature Tph in- 
stead of Tg. The quantity Tph is the electronic temper- 
ature commonly used in a self-consistent linear response 
calculation. Usually Tph is much larger then the phys- 
ical temperature Tq, (e.g. Tph « 4000 Kf"^. The main 
consequence of the variational property of the functional 
Fjj is that it allows us to compute phonon frequencies 
from an approximate force-constants matrix Cij{lj^Tq) 
in a non self- consistent way and performing an error that 
is quadratic in |n}(r, a;,To) — n}(r, 0, Tph)|. 

Thus, a time consuming self-consistent calculation of 
the non-adiabatic force-constants Ci.j{uj,Tq) has been 
replaced by an approximate non self-consistent one, 
Cij{uj, To), with an error that is negligible for the phonon 
frequencies, as it will be demonstrated in the applications 
considered in this work. 

We now first add and subtract C7j(0,rph) (the stan- 
dard linear-response adiabatic self-consistent force con- 
stants) in the left member of Eq. 21, and then perform 



where 



C,,(q,cj,To) = n,,(q,c^,To) 

+c,,(q,o,rph) 



n,,(q,c^,ro) = 



(22) 



E 



/k^(^o) - /k+qj(?o) 

df,(k,k + q)d!;,(k + q,k) 



Wfc(Tph) 



Affe(Tph) 



E 

kij 



.fkijTph) — /k+qj(^ph) 
f^ki ^k+qj 



d?-(k,k + q)d^,(k + q,k) 



(23) 



with the following definition of the deformation potential 
matrix element 



d™„(k + q,k) 



qrn\— \kn) 



Su, 



(24) 



qs 



and where |kn) is the periodic part of the Bloch wave- 
function, i.e. IV'kn) ~ e'^^''^\'kn)/^/Nk, and Uq^ is the 
Fourier transform of the phonon displacement u^,,,. The 
integration in Eq. 24 is understood to be on the unit 
cell. The quantity tiscF is the periodic part of the static 
self-consistent potential, namely VscF(r) = wscF(r)e*'i'". 
This static self-consistent potential is calculated using 
standard linear response at temperature Tph. 

The calculation of nsr(q, w, To) requires the knowledge 
of band energies and eigenfunctions on a denser Nk(To) 
k-point grid using an electronic temperature Tq only in 
an energy window of width ^ inax(rph,ci;) around the 
Fermi level. Moreover it docs not require to recalculate 
the derivative of the self-consistent potential. As such 
it is much less time consuming then an linear-response 
self-consistent calculation to obtain Crs(q, 0,ro)- The 
non-adiabatic phonon frequencies are calculated (in the 
clean limit) non self-consistently from the Hermitian part 

of C'sr(q, w. To) (see Eq. 5). 

The expression of the adiabatic force-constants at the 
physical temperature Tq is obtained setting w = in the 
Eq. 22: 



C,,(q,cj = 0,To) = n,,(q,0,ro) 
+C,,(q,0,Tph) 



(25) 



In this case the force constants are Hermitian and the 
error in the approximated force constants is quadratic in 
|n}(r,0,To)-n}(r,0,Tph)|. 

The advantage of the present procedure is that the 
linear-response self-consistent calculation is performed 
with a small number of k-points A^fe(Tph) whereas the 
low temperature force constants are obtained with a non 
self-consistent calculation over much denser NkiTo) k- 
points mesh. This approach requires, as in a conventional 



electron-phonon coupling calculation, the knowledge of 
the wavefunctions in a dense -/Vfc(To) k-points grid and 
in an energy window of the order of the maximum be- 
tween Tph and uj around the Fermi level, but does not 
require the self-consistent linear-response calculation of 
the derivative of the Kohn-Sham potential with respect 
to an atomic displacement. As such the procedure is 
substantially less time consuming. 



G. Phonon frequencies interpolation over k-points 

at fixed phonon-momentum: practical 

implementation 

The practical implementation of the above theoretical 
formulation for the phonon frequencies proceeds as fol- 
lows: 

1 . Perform a standard linear response calculation for 
a given phonon momentum q using a Nk{Tph) k- 
points mesh and smearing Tph for the electronic 
integration to obtain Csr(q, 0,Tph). 

2. In order to perform the second summation in 
nsr(q, i^,2o) (Eq. 23), calculate deformation po- 
tential matrix clement of Eq. 24 on the electron- 
momentum grid composed of A^fe(Tph) k-points us- 
ing a smearing Tph . 

3. Generate wavefunctions |kn) and energies etn on 
the denser iVfc(To) electron- momentum k-points 
grid. 

4. Perform a second non self-consistent calculation of 
the deformation-potential matrix-element on the 
more dense Nk{To) k-points grid using smearing Tq 
and wscF obtained on the Nk{Tp\i) grid to obtain 
the first summation in IlsrCq, ^^j To). 

5. Calculate C^rCq, w,ro) (or Csr{<l,^ = 0,ro)) using 
Eq. 22 (or Eq. 25). 

6. Diagonalize the dynamical matrix to obtain phonon 
frequencies and phonon eigenvectors. 

The procedure illustrated in this section is used to 
obtain at a fixed phonon momentum q well converged 
phonon frequencies with respect to electronic k-points 
and smearing. Then standard Fourier interpolation can 
be performed to obtain dynamical matrices throughout 
the BZ. However, this procedure, has still two main com- 
putational shortcomings. First, the calculations of the 
wavefunctions and matrix elements (points 3 and 4) be- 
come cumbersome when very dense NkiTo) k-points grid 
for the electron-momentum are necessary to converge. 
Second, in presence of phonon anomalies and not-too- 
smooth phonon dispersion a dense k-sampling of the 
phonon BZ is required. 

An optimal solution to overcome these problems is 
represented by the use of maximally localized Wannier 
functions^^'^* as an alternative electron basis function. 



The method and the implementation of the Wannier 
functions approach for the electron-phonon interpolation 
is presented in Sec. III. 



H. Phonon-linewidth and of the electron-phonon 
coupling 

The imaginary part of the force-constants matrix is re- 
lated to phonon-damping, namely the energy-conserving 
decay of a phonon in particle-hole pairs. 

The imaginary part of Eq. 16 is: 
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where P indicates the principal part. To obtain Eq. 26 
we have assume that the unperturbed Hamiltonian is 
time-reversal symmetric so that a real set of eigenfunc- 
tions exists. 

This expression is exact but it has two disadvantages. 
First it requires the expensive calculation of the deriva- 
tive of the self-consistent potential at finite frequency 
both in its real and imaginary parts. Second the phonon- 
decay in electron-hole pairs is not manifest in the term 
including the principal part. These problems persist if 
the imaginary part of Eq. 17 is considered, since, in 
this functional, both (5VscF(r,a;)/(5u/ and n}(r, li;,T) are 
complex quantities. These two shortcomings are absent 
if the approximated force-constants matrix in Eq. 21 is 
used. Indeed the imaginary part of Eq. 21 is: 
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with an error quadratic in \n\{r,ijj) — n}(r, 0)|. Fourier 
transforming and replacing Eq. 27 in Eq. 9 gives the well 
known equation for the phonon linewidth 7q^ full-width 
half- maximum (FWHM), that is: 
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In absence of the electron-electron interaction a similar 
expression involving only the derivative of the bare exter- 
nal potential is easily obtained using Fermi golden rule. 
In the presence of electron-electron interaction our for- 
malism justify the replacement of the derivative of the 
bare external potential by the derivative of the static 
screened potential VscF(r,0) in the Fermi golden rule, 
with a well controlled error. 

Under certain conditions illustrated in Refs. 29,30, the 
phonon linewidth can be reduced to: 
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Finally, neglecting ujq^ in Eq. 29, we obtain Allen 
formula^^ namely, 
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(30) 
The AUen-formula^'^^ is widely used, and represents a 
good estimation of the phonon-linewidth due to electron- 
phonon effects, in the absence of anharmonic effects and 
other scattering processes. In addition, the Allen-formula 
relates the phonon-linewidth, as measured by inelastic 
Neutron or X-ray measurements, to the electron-phonon 
coupling as: 






(31) 



where A/'s is the electronic density of states per spin at 
the Fermi energy. 

With this definition, the isotropic Eliashberg-function 
is 
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where Nq is the number of phonon wavevectors. We also 
define the integrated electron-phonon coupling as 
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III. WANNIER FUNCTIONS 



Wannier interpolation of the electron-phonon 
matrix element 



In this section we explain how to interpolate the 
electron-phonon matrix element throughout the BZ using 
Wannier functions^". 

Using standard first-principles methods a set of Bloch 
functions V'kn are generated with k belonging to a uni- 
form N^ k-points grid centered in T. A set of Wannier 



functions centered on site R are defined by the relation 
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A suitable transformation matrix, t/m„(k), must be de- 
termined in the so-called Wannierization procedure. In 
this work we choose to work with Maximally Localized 
Wannier functions (MEWF)^**'^^, although other Wan- 
nierization schemes are possible. In the MLWF case, the 
matrices C/„m(k) are obtained following the prescription 
of Refs. 14,28, which guarantees the maximal space lo- 
calization of the final Wannier functions. This last re- 
quirement will be essential in the spirit of interpolation 
of the electron-phonon matrix elements. 

Specifically, if the band-structure is formed by a com- 
posite group of bands and the number of Wannier func- 
tions is identical to the number of bands of the composite 
group of bands, then the matrix L'^mn(k) is a square ma- 
trix. On the contrary if the desired bands are entangled 
in a larger manifold, a preliminary disentanglement pro- 
cedure from the manifold must be performed and then 
the square matrix f7mn(k) is obtained^^. 

The deformation potential matrix-elements d^„(k -|- 
q, k) (see Eq. 24) are calculated using linear response 
scheme. In this calculation particular care is needed 
for the deformation potential at zone center, namely 
d^„(k, k), as explained in sec. IIIB. 

It is crucial to note that the periodic parts jkn) and 
jk + qm) entering in d^„(k-|-q, k) have to be exactly the 
same wavefunctions used for the Wannierization proce- 
dure. If this is not the case, spurious (unphysical) phases 
appear in the |kn)'s (i.e. a phase added by the diago- 
nalization routine or other computational reasons) and 
the localization properties of the Wannier functions in 
real space are completely lost. A way to enforce this 
condition is to use a uniform grid centered at F so that 
k + q = k' + G with k' still belonging to the original 
grid. In this way jk+qn) can be obtained from |k'n) sim- 
ply multiplying by a phase determined by the G-vcctor 
translation and no additional phases occur. 

Exploiting translational invariance, the deformation- 
potential matrix-element in the Wannier function basis 
is obtained by Fourier transform as. 
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where R and R^ belong to a TV™ real-space supercell. 
At this point it is important to underline that the grid 
of N^ k— points on which the Wannierization process 
has been carried out has to be exactly the same grid for 
the phonon-momentum on which the dynamical matrices 



have been computed. One linear response calculation for 
each k-point in the phonon irreducible BZ (IBZ) of the 
NJf electron k-points grid has to be carried out. If this 
constraint is not respected, the localization properties of 
the electron-phonon matrix element in real space is not 
guaranteed. 

Inverting Eq. 34, we obtain 
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Noting that d^„(k+q, k) = ^(^k+qml^l^lV'kn), and 

k *-l® 

using Eq. 36 one obtains: 
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where R and R^ belong to a iV^ real-space supercell. 
Now, if d^„(R, Rl) is localized inside the NJ^ real-space 
supercell then the interaction between different N]f real- 
space supercells can be neglected and df^^Ck + q, k) can 
be obtained from Eq. 37 via a slow Fourier transform. In 
practice this means that in Eq. 37 now k and q are any k- 
points in the BZ. As a result of the interpolation, d^„(k+ 
q, k) can be used to calculate any physical property as 
in a simple tight-binding scheme. 



is well defined and independent on the direction of q. In 
general one has Aq+s y^ 0. 

At q = 0, the calculation is performed in the canonical 
ensemble, with a constant number of electrons. In this 
calculation the variation of the average Coulomb poten- 
tial is conventionally set to zero, namely: 
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To keep constant the number of electrons the deriva- 
tive of the Fermi energy with respect to the phonon 
displacement SeF/SuOs can be different from zero. In 
particular^'^^ it must be: 
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In the QUANTUM-ESPRESSO implementation^! the 
<5i^scF(r)/(5uqs stored-on-disk is discontinuous at q = 0, 
since, in this case, it does not include (as it should) the 
contribution from the average Coulomb potential. As a 
consequence the electron phonon matrix elements com- 
puted with (5wscF(r)/(5uos arc incorrect. Moreover such 
discontinuity deteriorates the localization properties of 
the electron-phonon matrix elements in real space. The 
problem is easily solved by redefining the self-consistent 
potential to be used in the calculation of the deformation- 
potential matrix-elements for q = as: 



B. Deformation potential matrix elements of 
optical zone-center phonons 
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In this section we discuss the peculiarities related to 
the calculation of the electron-phonon matrix elements 
for optical zone center phonons (q = 0). The peri- 
odic part of the self-consistent potential, induced by 
the phonon displacement at wavelength q of the atom 
s can be decomposed in the Coulomb (CI) potential (the 
sum of the bare and Hartree potentials) and exchange- 
correlation (XC) contributions: 
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The integral of the Coulomb contribution over the unit 
cell volume fl is: 
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In DFT linear response codes, such as QUANTUM- 
ESPRESSO^!, the phonon calculation with q 7^ and 
q = are treated with two different approaches. ^'^^ 

At q ^ 0, the calculation is performed within the gran- 
canonical ensemble, with a constant electron chemical- 
potential (Fermi level) ep. In a metal, the limit 



lim Aqs = Aq+s 



(40) 



Note that since Sep/Suos transforms under symmetry op- 
eration as a force, it is different from zero only for the 
atoms for which the internal coordinates are not fixed by 
symmetry. For this reason, Sep/Suos is finite in CaCg, 
while it vanishes in MgB2. 

We remark that the same problem occurs when a 
frozcn-phonon calculation of ^|^''^'"^ is carried out as in 
all the electronic-structure codes the average Coulomb 
potential is conventionally set to zero. 



C. Wannier interpolation of the dynamical matrix 
at fixed phonon-momentum. 

Steps 4, 5, 6 of section II G can be repeated using Wan- 
nier functions to interpolate at fixed phonon momentum 
q the quantity Ilsr{(l,Lo,To) to obtain Eq. 23. The ad- 
vantage of the Wannier- interpolation-based strategy with 
respect to the procedure explained in sec. II G is that the 
time needed to calculate the matrix element is now neg- 
ligible as there is no need to obtain, from first-principles, 
wavefunctions and energies on the dense Nk{To) k-point 
grid. Still it is necessary to perform Fourier interpolation 
to obtain dynamical matrices throughout the full BZ. 



D. Wannier interpolation of the dynamical matrix. 



In what follows we outline a strategy to obtain dynam- 
ical matrices throughout the BZ using both Fourier inter- 
polation and Wannier functions. The method solves the 
main problem related to the presence of Kohn-anomalies, 
originating from long range interactions in the dynamical 
matrix. The idea is to separate in the force-constant ma- 
trix obtained in Eq. 22 in the short and long range com- 
ponents. The long range force constants are associated 
to Kohn anomalies driven by Fermi surface nesting, and, 
as such, they cannot be easily Fourier interpolated and 
require an accurate sampling of the Fermi surface. This 
last part will be treated using the Wannier interpolation 
scheme explained in section III A. On the contrary short 
range force constants produce smooth phonon dispersions 
and can be safely Fourier interpolated and do not require 
accurate sampling of the Fermi surface. Treating sepa- 
rately the two length scales permits to obtain converged 
phonon frequencies with respect to k-point sampling and 
electronic temperature, everywhere in the Brillouin zone. 

This procedure is similar to the treatment of long- 
range dipolar forces due to the Born effective charges 
in polar semiconductors^'^. In this case, these long range 
forces generate a non-analytical behavior of the dynam- 
ical matrix at zone center. While Fourier interpolation 
of the dynamical matrices including the non-analytical 
terms is impossible, it becomes possible if the long- 
range forces are subtracted to obtain short-range force- 
constants. The non-analytical part subtracted to per- 
form Fourier interpolation is added again after Fourier 
interpolation. 

In this spirit, we rewrite Eq. 22 as, 

a,(q,w,ro)=n,,(q,c^,To) 
-n,,(q,0,Too) + C,,(q,0,T„o) (44) 

where 

a,(q,o,roo) = C,,(q,0,rph) + 

+ n,,(q,0,T„o) (45) 

and Too is an electronic temperature large enough in 
order to have only short range force constants named 
Csr(q, 0, Too), and no Kohn anomalies in the correspond- 
ing phonon branches. As a consequence, Fourier interpo- 
lation can be applied to Csr i<l, 0, Too). 

In practical implementation the procedure is the fol- 
lowing: 

1. Obtain Csr(q, 0,Tph) from self-consistent linear- 
response phonon calculations. Such force constants 
are evaluated using a N]f k-point grid for the 
phonon momentum and an iVfe(Tph) k-point grid 
for the electron momentum. 

2. Calculate ns-r(q, 0,Too) using the Wannier func- 
tions as described in section IIIC. 



3. Use Eq. 45 to compute C'sr(q,0,Too) on the N^ 
phonon momentum grid. 

4. Fourier interpolate Csr(q, 0,Too) to obtain dynam- 
ical matrices at any desired phonon momentum at 
temperature Too- 

5. Obtain C'sr(q, w,To) in Eq.44 by Wannier interpo- 
lation of nsr(q,w. To)— nsr(q, 0, Too) at any desired 
phonon momentum. 

6. Diagonalizc the dynamical matrix associated to 
C'sr(q, w,To,Tph) to obtain adiabatic, w = 0, or 
non-adiabatic (in the clean limit) w = Wq^, phonon 
frequencies. These phonon frequencies include 
Kohn-anomalies as long range terms are properly 
treated by the Wannier-functions approach. 



IV. APPLICATIONS 

We applied the theoretical and computational frame- 
work developed in the preceding sections to two well 
known superconductors, namely magnesium diboride 
(MgB2) and Calcium intercalated graphite (CaCg). In 
these systems phonon, electron-phonon and supercon- 
ducting properties are extensively investigated from both 
theoretical and experimental point of view and a detailed 
comparison can me made. At the same time, they rep- 
resent non-trivial examples, where convergence problems 
in first-principles approaches are relevant for the calcu- 
lation of phonon and superconducting properties. 



A. MgB2 

Magnesium diboride (MgB2) is, without doubts, one 
of the most studied and investigated materials, since 
2001. The discovery of superconductivity at interme- 
diate temperatures (40 K)^^, in this compound, repre- 
sented an unexpected surprise in the scientific commu- 
nity. Easily available, a huge amount of experimental 
characterizations were immediately possible^'^. On the 
other hand, MgB2 being composed of light elements in 
an hexagonal unit cell with only three atoms, allowed 
a first-principles computational approach for the theo- 
retical description and understanding of its normal and 
superconducting properties. Density Functional Theory 
and linear response techniques, revealed an extraordinary 
high electron-phonon coupling between sp'^ bonding and 
anti-bonding boron orbitals, with a particular in-plane 
mode {E2g) of the hexagonal boron sheet'^^"^^. Based 
on model calculations of the electron-phonon coupling 
parameter (A ~ 1 in Ref. 35), the approximate criti- 
cal temperature was estimated to be around 40 K. How- 
ever, a subsequent careful investigation of the problem, 
revealed that the first-principles calculated values of A, 
were too low to explain this high critical temperature, 
unless unphysical low Coulomb pseudo-potentials were 
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used'^^. The solution of this inconsistency, has come with 
the discovery of two-band-superconductivity'^* in MgB2. 
Two different gaps develop on the a and tt bands of the 
boron sheets'^^. This allows to gain new scattering chan- 
nels to enhance Tc still containing the average value of 
A. At the moment, MgB2 is classified as a two-band 
electron-phonon superconductor'^^^^^ . 

This continuous improvement of the theoretical expla- 
nation of the normal and superconducting properties of 
MgB2 was in part due to some difficulties in the proper 
calculations of the electron-phonon coupling and of the 
phonon frequencies. Soon after the discovery of super- 
conductivity in MgB2, the first theoretical papers (based 
on standard DFT codes) evidenced a relevant and anoma- 
lous difficulty to properly "converge" the Y,2g phonon 
frequency and the electron-phonon coupling. The rea- 
son is now clear: due to the two dimensional nature 
(almost cylindrical) of the a* Fermi surfaces electron- 
phonon coupling calculations require a careful integration 
of the Fermi surface is needed. This translates, in prac- 
tical calculations, with the necessity to use an extremely 
dense k- space integration for electronic and phononic 
degree of freedom''^. Although feasible in principle, it 
would require a huge amount of resources and computa- 
tional time. 

So, still after about ten years of calculations, the funda- 
mental question about the value of the electron-phonon 
coupling in MgB2, remains unanswered. This informa- 
tion is fundamental in order to validate the electron- 
phonon origin of the superconductivity in MgB2 and at 
the same time, include other, still unexplored pairing in- 
teractions. 

The Wannnier interpolation scheme, developed in 
the present paper, of phonon frequencies and electron- 
phonon matrix elements, represents the main tool to 
solve this problem. It retains the first-principles char- 
acter, but at the same time it is affordable from a com- 
putational point of view. 

In MgB2 the BZ-center phonon-frequency, measured 
by Raman spectroscopy, lies between the adiabatic and 
clean-limit-non-adiabatic phonon-frequency^^ computed 
by DFT. Thus MgB2 is in the intermediate regime (see 
Sec. II B and the relaxation time estimated in Ref. 12). 
Since the experimental frequency at BZ center is closer to 
the adiabatic than the clean-limit-non-adiabatic value, ^^ 
in MgB2 we will not consider non-adiabatic effects, and 
we will present only the adiabatic dispersion. 



1. Technical details 

Wannier-interpolation of the MgB2 band-structure is 
carried out using a NJ^ = 6x6x4 k-points grid centered 
at r using 5 Wannier functions starting from random 
projections. The Wannierized band structure is in excel- 
lent agreement with first-principles calculation in a 1.5 
eV region from the Fermi level, as shown in Fig. 1. 

The phonon dispersion and the electron-phonon ma- 




FIG. 1: (color online) Wannier-interpolated (red-dashed) and 
first-principles-calculated bands (black-continuous) for MgB2 . 
The band-energies are plotted with respect to the Fermi level. 



trix element are calculated using linear-response ' on 
the irreducible part of the NJ^ k-points grid centered at 
F. The 5vscF/Suqs on the full grid are obtained ap- 
plying symmetry operations of the crystal. In the linear 
response calculation we perform electronic integration us- 
ing a Nk = 16x16x12 fc— point grid in the Brillouin zone 
and an Hcrmite-Gaussian smearing Tph — 0.025Ryd. 

The interpolation of phonon frequencies is performed 
along the strategy outlined in sec. Ill D and using Too — 
0.11 Ryd to obtain smoothed "high temperature" phonon 
frequencies. 



2. Adiabatic phonon dispersion 

The Wannier-interpolated adiabatic phonon dis- 
persions on different k-point grids and smearings 
are compared to state-of-the-art electronic structure 
calculations'^^'^^''''^ and experimental data^'*^ in Fig. 2. 
In the top panel of Fig. 2 we compare the Wannier in- 
terpolated phonon dispersion with Tq = 0.02 Ryd and 
Nk{To) = 30'^ with linear response calculations per- 
formed on selected points in the Brillouin zone using the 
same mesh and the same electronic temperature. The 
error in the interpolation scheme is always smaller then 
0.5 meV. In the bottom panel of Fig. 2 we compare Wan- 
nier interpolated dispersion curves on Nfe(To) = 45'^ at 
To — 0.01 Ryd with Fourier interpolated linear response 
calculations from a 6 x 6 x 4 phonon momentum grid. As 
can be seen anomalies are clearly found along TA using 
our new method but are completely missed by standard 
Fourier interpolation. The main differences with respect 
to previous calculations'^^'^^'"''^ are : 

1. a prominent softening (Kohn anomaly) of the two 
E2g modes at q ~ 0.25rK is seen. The Kohn 
anomaly is also present along the FM high sym- 
metry direction and in a circular region around 
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FIG. 2: (color online) Top panel: Wannier-interpolated adi- 
abatic phonon dispersion (red line) compared to standard 
linear-response calculations performed on selected points in 
the Brillouin zone (black dots). Bottom-panel: Comparison 
between state-of-the-art calculated phonon dispersion using 
Fourier interpolation (black line) and Wannier interpolation 
(red line). Both interpolation schemes start from the same 
linear response calculation. Experimental data in the bot- 
tom panel are from Ref. 5 (blue circles) and from Ref. 44 
(black circles). Electronic temperatures are in Rydberg. The 
acoustic sum rule is not applied. Interpolated branches are 
obtained by connecting the points at which the interpolation 
has been carried out in the order of increasing energy. 



the zone center. Its origin is due to a in-plane 
nesting between the two a cylinders. This Kohn 
anomaly could be probably inferred from previ- 
ous calculations^'''^'^'^ on smaller k-points grids al- 
though the softening of the E2g modes in these 
works is weaker. 

2. A second Kohn anomaly, absent in previous 
calculations^''*^''*'^, is present on the E2g and Big 
modes along the kz direction at q w O.SFA. This 
Kohn anomaly is due to nesting between differ- 
ent cylinders along FA. As it can be seen by the 



comparison with previously available experimental 
data in Fig. 2, its presence is confirmed by inelas- 
tic X-ray scattering measurements of the phonon 
dispersion^'^''. The data at q « O.SFA are indeed 
in disagreement with previous phonon dispersion 
calculation^'^^"'*^"'*"' even if this disagreement was 
overlooked in all previous publications. 

3. Overall the Wannier-interpolated phonon disper- 
sion is in better agreement against experiments 
with respect to the linear-response calculated one. 

These differences point out the need of having an ultra- 
dense k-point sampling of the Brillouin zone not only 
for the electron-phonon coupling but also for the phonon 
dispersion. 



3. Electron-phonon coupling 

Having interpolated the electron-phonon matrix ele- 
ments and the phonon frequencies it is possible to calcu- 
late the phonon-linewidth and the electron-phonon cou- 
pling with a high degree of precision, eliminating source 
of errors coming from insufficient fc— point sampling. In 
Fig. 3 the Wannier-interpolated Eliashberg-function is 
compared with previous calculations done with different 
approaches. There are two main improvements due to 
the better k-point sampling. 

First, several previous calculations^^'^^'"'^ on smaller 
k-point grids generally overestimate the position of the 
main peak in a^Fiuj) (centered around cj « 63 meV 
in our work) due to the coupling of cr-bands with the 
E2g phonon modes. This happens because in previous 
calculations the Kohn anomalies of the E2g branches 
at q w 0.25Fiir was underestimated and because the 
anomaly along TA is missing. This crucial dependence 
of the Eliashberg function on the E2g phonon frequency 
demonstrates the need of having both an accurate deter- 
mination of A and of the phonon frequencies to describe 
superconducting properties. 

Second, in our work there is substantial more weight 
then in previous works^^'**^ in the 90-100 meV region. 
This energy region is dominated by the coupling of tt 
electrons to the E2g vibrations. 

We now compare in more details our work with previ- 
ous calculations. In the work of Bohnen et al.^^ all the 
peaks are up-shifted not only by the limited k-points sam- 
pling but also by the choice of using the DFT-minimized 
crystal structure and not the experimental one. Kong 
and coworkers'*^ have calculated \qy and Wq^ on a 6 x 6 x 6 
phonon-momentum grid. Then for q of this grid be- 
longing to the FA direction the corresponding Aqjy val- 
ues were replaced with some of the 12 x 12 x 6 k-point 
grid (see Ref. 42 for more details). Fig. 3 shows that 
this approximation strongly overestimates the energy of 
the peak due to coupling between a electronic states and 
the E2g mode. Indeed the average electron-phonon cou- 
pling is much larger than the majority of the calcula- 
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TABLE I: Comparison between different electron-phonon cou- 
phng calculations in MgB2 using different first principles 
methods. The label Nfc (Ng) indicates the number of k-point 
used for electron (phonon) momentum integration. 



the phonon dispersion does not present the two observed 
Kohn anomalies and the main peak of the a'^F{ijj) is at 
higher energies respect to our work. This difference re- 
flects the need of having an accurate determination of the 
phonon frequency in order to have converged Ehashberg 
functions. This appears in the lower phonon frequency 
logarithmic average obtained in our work with respect to 
previous calculations as shown in table L 



FIG. 3: (color online) Comparison of our Wannier- 
interpolated MgB2 Eliashberg function with previous calcu- 
lations. Kong et al.'^'^ used first-principles calculations and 
a particular treatment for Aq^g^ for q along FA. Bohnen et 
al.*^ used linear response calculations while Choi et al.'^^ per- 
formed frozen-phonon calculations. Finally Eiguren et al.^^ 
used a different implementation of Wannier interpolation on 
the electron-phonon matrix elements with 40'^ k-points grid 
for electron-momentum and 40^ k-points grid for phonon mo- 
mentum, but did not interpolate phonon frequencies. In this 
work we Wannier-interpolate both electron-phonon coupling 
(using 80'^ k-points grid for electron-momentum and 40^ k- 
points grid for phonon momentum) and phonon frequencies 
(using 30^ k-points grid for electron-momentum and 40^ k- 
point grid for phonon momentum). The integrated A(a;) from 
our work is also shown in the bottom panel. 



tions present in literature (see Table I). Choi et al.^^ 
used the LDA-minimized lattice structure but surpris- 
ingly obtained softer E2g phonon frequencies at F then 
all other works using experimental lattice parameters. In 
Ref. 39 the electron-phonon coupling is calculated using 
an interpolation scheme that does not assure the local- 
ization properties of the electron-phonon matrix element 
in real space. Nevertheless the position of the main peak 
in a'^F{uj) is similar to what we obtain. On the con- 
trary the secondary peak related to coupling between tt 
states and the E2g mode is at too low energy. Finally, 
Eiguren et al.^^ used a Wannier interpolation scheme 
that requires an explicit calculation of the wavefunctions 
for any phonon momentum. The interpolation scheme 
was used to calculate the average electron-phonon cou- 
pling but not to interpolate phonon frequencies as they 
were obtained by Fourier interpolation. As a consequence 



B. CaCfi 

Graphite intercalated compounds represent a wide 
class of compounds. Due to the particular bonding prop- 
erties of carbon atoms, with strong and stable sp^ bonds 
and weak Van der Walls tt^ bonds, graphite allows incor- 
poration of many different atomic species between the 
hexagonal C sheets. Alkali metals and alkaline-earth 
metals are the most common, but even H, Hg, Tl, Bi, 
As, O, S and halogens are reported^^ 

The possibility to change the number of electrons in 
the covalent C-C bonds, always stimulated the search of 
superconductivity in intercalated graphite, but only in 
2005 calcium intercalated graphite (CaCe) was found to 
superconduct at temperatures as high as 11.5 K'^^. 

First-principles theoretical investigations on CaCe re- 
vealed that the superconducting critical temperature 
can be explained considering only the electron-phonon 
mechanism^^'^* , even if other pairing mechanisms were 
proposed (excitonic and plasmonic). 

The pairing originates from the electron-phonon in- 
teraction between intercalant Fermi Surface and the C 
out-of-plane modes and Ca in-plane modes*^. These last 
modes, with very low frequencies (~ 15 meV), account 
for about one- half of the total A. 

In this section, we re-examine the phonon and electron- 
phonon properties of CaCg, by means of the computa- 
tional framework based on the Wannier interpolation. 

In fact, from a computational point of view, the cal- 
culation of the energy position and q-dispersion of the 
low energy Ca modes requires very high precision, not 
obtainable with the state-of-the-art phonon and electron- 
phonon calculations. 
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FIG. 4: (color online) Wannierized (red-dashed) and first- 
principles-calculated bands (black-continuous) for CaCe- The 
horizontal dashed line marks the Fermi level. 



CaCg is well described by the clean limit regime. In- 
deed, the CaCe zone-center phonon frequency measured 
by Raman, -'^^'''^'^° is significantly larger than the DFT 
adiabatic one but agrees very well with the DFT clean- 
limit- non-adiabatic frequency^^'^^'"'^''^". For this reason 
we will discuss both the adiabatic and the non-adiabatic 
dispersion. 



Technical details 



Wannier-interpolation of the CaCg band-structure is 
carried out using a A^^ = 6x6x6 k-points grid centered 
at r. We use 7 Wannier functions starting from an initial 
configuration with one s-state on Calcium and one Pz 
state on each Carbon. The Wannierized band structure 
as compared to first-principles calculations is shown in 
Fig. 4. 

We perform one linear response calculation for 
each point of the irreducible N^ A:— point grid using 
^fe(Tph) = 8x8x8 and Hermitian-Gaussian smearing 
Tph = 0.05 Ryd. We then obtain Svscp/Suqs on the full 
grid applying the symmetry operations of the crystal. At 
zone center, we added a Fermi level shift to the q = 
component of SvscF/Suqs to obtain an analytical behav- 
ior of at q = (see section IIIB for more details). 

The interpolation of phonon frequencies is performed 
along the strategy outlined in sec. HID and using 
Too — 0.075 Ryd to obtain smoothed "high temperature" 
phonon frequencies. For the calculation of non-adiabatic 
phonon frequencies we chose w in Eq. 21 as the E2g 
phonon frequency at F and we also use A^fc(To) = 50^^, 
To = 0.01 Ryd and r? = 0.01 Ryd in Eq. 23. 



2. Adiabatic phonon dispersion 

The phonon dispersions of CaCg as obtained from 
Fourier-interpolation and from Wannier-interpolation are 
illustrated in Fig. 5. The Fourier interpolate phonon 
dispersion is in agreement with previous linear re- 
sponse calculations^^'^^'^^. When compared with the 
more converged calculations performed using Wannier- 
interpolation several differences are found. The largest 
deviation concerns the high-energy modes at sa 160— 185 
meV, in particular along TL. These modes involve car- 
bon vibrations parallel to the Graphite layers {Cxy vibra- 
tions). There is a substantial hardening of the two high- 
est optical phonon-modcs approaching zone-boundary 
and a somewhat weaker softening close to zone center. 
The shape of the phonon dispersion of the two high- 
est optical modes along TL recalls the typical behav- 
ior of the real-part of the phonon self-energy due to 
electron-phonon coupling (see for example for the case 
of MgB2 Ref. 30, fig 3a.). The electron-phonon coupling 
to these optical modes comes essentially from tt* states 
and mostly from intraband-coupling'^'^ (small Fermi sur- 
face) so that a quite accurate BZ sampling is necessary to 
display this behavior. This explains the strong fc— point 
dependence of the result. 

Moreover the relevance of having an accurate sampling 
of the Fermi surface is evident from the large number of 
Kohn-anomalies (arrows in Fig. 5 mark the occurrence 
of a Kohn-anomaly) present in the Wannier-interpolated 
branches and almost always absent in the Fourier inter- 
polated ones. At high energy (150-180 meV) we count 
a large number of anomalies, no one of which is present 
in Fourier interpolated branches. These anomalies are 
mostly located close to the X and x points. At low en- 
ergy we count 8 more softenings, 6 of which involve points 
close to the X-point (although on different branches) and 
two occur at the x point. In Fourier interpolated calcu- 
lations only a very broad anomaly is visible in the low 
energy modes, located exactly at the X-point, nothing 
in other energy modes. The facts that close to the X 
point (i) several anomalies at the same vector occur on 
different modes and that (ii) the phonon frequencies as- 
sociated to the anomaly are strongly dependent on the 
electronic temperature suggest an origin dependent from 
the Fermi surface geometry. Indeed this is the case and 
the anomaly is illustrated in the nesting vector in Fig. 6. 



3. Non-adiabatic phonon dispersion 

Non-adiabatic effects occur in CaCg close to zone 
center^^'"'^^. However it is unknown how strong are these 
effects far from zone center, namely what is the degree of 
localization of non-adiabatic effects around the F point. 
The question is meaningful as if non-adiabatic effects ex- 
tend substantially in the Brillouin zone then they can 
be relevant for superconductivity and transport as the 
average electron-phonon coupling and the critical tem- 
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FIG. 5: (color online) Wannier-interpolated (red-continuous) and Fourier-interpolated (black-dashed) adiabatic phonon dis- 
persion in CaCe. The center and right panels are zooms to the high and low energy parts of the phonon dispersion. The 
Fourier-interpolated phonon branches are obtained from the dynamical matrices calculated on a NJf = 6^ phonon grid, using 
Nk = 8 and Hermite-Gaussian smearing Tph = 0.05 Ryd. The Wannier interpolated bands are obtained on a A^^ — 40 and 
Hermitian-Gaussian smearing Tph = 0.01 Ryd. The Kohn-anomalies are indicated by arrows. 



Reference 


A 


a;iog(meV) 


Nk 


N, 


Calandra et al. [47] 


0.83 


24.7 


6x6x6 


4x4x4 


Kim et al. [52] 


0.84 


26.28 


8x8x8 


4x4x4 


Sauna et al. [48] 


0.85 


28.0 


6x6x6 


6x6x6 


This work 


0.829 


27.9 


59^ 


20^ 



TABLE II: Comparison between different electron-phonon 
coupling calculations for CaCe using different first principles 
methods. The label N^ (N,) indicates the number of k-point 
used for electron (phonon) momentum integration. 



perature could be affected. 

In Fig. 7 we plot the Wannier-interpolated non- 
adiabatic phonon dispersion. At F we obtain w'^^ ~ 190 
meV (« 1532 cm"^), in good agreement with what found 
in Ref. 12. Most interestingly the full E2g branch is 
shifted to higher energies (particularly along FL, the k^ 
direction) by the non-adiabatic effects. Even at half-way 
to the zone border the correction is sizable and it can be 
easily measured with inelastic X-ray or Neutron scatter- 
ing experiments. In the case of CaCe as the high-energy 
modes are weakly-coupled to electrons the non-adiabatic 
correction, even if large, probably does not affect the 
critical temperature. The situation can be substantially 
different in other superconductors so that non-adiabatic 
effects have to be taken into account to describe the 
electron-phonon interaction in layered superconductors. 



4- Electron-phonon coupling 

The Eliashberg function calculated with the adiabatic 
phonon frequencies and the integrated electron-phonon 
coupling for CaCg are illustrated in Fig. 8. The situation 



is very similar to previous works in what concern the 
average value of the electron-phonon coupling, as shown 
in table II. 

Previous calculations performed using standard linear- 
response methods report similar values of A and a 13% 
variation is present on wiog (see table II). Our calculation 
confirms the largest value reported in literature. 



V. CONCLUSIONS 

In this work we developed a first-principles variational 
scheme to calculate adiabatic and non-adiabatic phonon 
frequencies in the full Brillouin zone. We demonstrated 
how the method permits the calculation of phonon dis- 
persion curves free from convergence issues related to 
Brillouin zone sampling. Our approach also justify the 
use of the static screened potential in the calculation of 
the phonon lincwidth due to decay in electron-hole pairs. 

We apply the method to the calculation of the phonon 
dispersion and electron-phonon coupling in MgB2 and 
CaCg. In both compounds we demonstrate the oc- 
currence of several Kohn anomalies, absent in previous 
works, that are manifest only after careful electron and 
phonon momentum integration. In MgB2, the presence of 
Kohn anomalies on the E2g branches improves the agree- 
ment with measured phonon spectra and affects the po- 
sition of the main peak in the Eliashberg function. In 
CaCg we show that the non-adiabatic effects on in-plane 
carbon vibrations are not localized at zone center but 
are sizable throughout the full Brillouin zone. Thus NA 
effects can be relevant for the determination of the su- 
perconducting properties of layered superconductors. In 
general our work underlines the important of obtaining 
accurate (well converged) phonon frequencies to deter- 
mine the electron-phonon coupling and superconducting 
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FIG. 6: (Color online) Fermi surface of CaCe calculated using Wannier-functions on a 100"^ k-points grid. The green arrows 
mark the nesting vectors responsible for the Kohn-anomaly close to X. 



properties in a reliable way. 

The large gain in the computational time neces- 
sary to obtain phonon dispersion provided by our the- 
oretical scheme opens new perspectives in large-scale 
first-principles calculations of dynamical properties and 
electron-phonon interaction. Furthermore, the varia- 
tional property of the force constants functional with 
respect to the first-order perturbation of the electronic 
charge density is a property that can be generalized to 
the calculation of other response functions such as op- 
tical properties and the calculation of Knight shifts as 
measured in nuclear magnetic resonance^^. 
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